Translation invariant tensor product states in a finite lattice system 
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We show that the matrix (or more generally tensor) product states in a finite translation invariant 
system can be accurately constructed from the same set of local matrices (or tensors) that are 
determined from an infinite lattice system in one or higher dimensions. This provides an efficient 
approach for studying translation invariant tensor product states in finite lattice systems. Two 
methods are introduced to determine these size-independent local tensors. 
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I. INTRODUCTION 

Periodic boundary conditions are useful for simulat- 
ing a large system by modeling a small part that is far 
away from its edges. In particular, a periodic system is 
free of boundary effects. It is easier to carry out finite 
size scaling analysis from the results obtained on peri- 
odic systems. This makes the extrapolation to the ther- 
modynamic limit more transparent and smaller systems 
are needed in the simulations. Moreover, the energy- 
momentum dispersion of excitation states can be better 
studied in translation invariant periodic systems. 

However, many numerical renormalization group 
methods, such as the density-matrix renormalization 
gropu (DMRG) 1 , works better in systems with open 
boundary conditions than those with periodic bound- 
ary conditions. This is because the DMRG wavefunction 
does not have the right entanglement structure in a pe- 
riodic system. The criterion for the basis truncation in 
the DMRG is governed by the the bipartite entanglement 
entropy, which is bounded by the logarithm of the basis 
number D retained in the truncation^. To calculate the 
bipartite entanglement entropy, one needs to split the 
system into two parts by cutting one bond in an open 
chain, but two bonds in a periodic chain. Thus the en- 
tanglement entropy in a one-dimensional (ID) periodic 
system is about twice that in the corresponding open 
system. In the DMRG calculation with open boundary 
condition, the computational cost scales as 0(LD 3 ) for a 
system of size L. However, to achieve the same accuracy, 
a periodic system needs roughly D 2 states per block. The 
computational cost scales roughly as 0(LD e ). In two or 
higher dimensions, the entanglement entropy grows faster 
with the system size. The number of states needed scales 
exponentially with the number of boundary sites if peri- 
odic boundary condition is imposed. 

The DMRG can be rephrased as a variational method 
over a class of matrix product states^ The matrix prod- 
uct states and related algorithms have been actively ex- 
plored. In 2004, Verstraete et al. pointed out that the 
matrix product wavefunction can also be used to cure 
the problem met in the DMRG with periodic boundary 
condition. 4 They proposed a variational approach to eval- 
uate the matrix product state and showed that the local 
matrix elements can be determined by solving a general- 
ized eigenvalue problem of dimension D 2 . The computa- 



tional cost of their algorithm scales as 0(LD 5 ). 

Recently, Pippan et al. proposed an approximate 
scheme to evaluate the variational matrices used in the 
generalized eigenvalue equation in terms of a singular 
value decomposition-. Their scheme reduces the compu- 
tational effort from 0(LD 5 ) to 0{LD 3 ), which is compa- 
rable to that of DMRG with open boundary condition. 
More recently, Pirvu et al. showed that the efficiency for 
evaluating matrix product states with periodic bound- 
ary condition can be further improved if the system is 
translation invariant £ The cost of this method scales as 
0(mD 3 ), where to is a number much smaller than L 
and becomes a constant for very large system size. A 
0(D 3 ) cost can be also achieved by applying the varia- 
tional Monte Carlo to a matrix product stated 

In this work, we further explore physical properties of 
translation invariant matrix or tensor product states in 
one or higher dimensions. We will show that one can use 
a set of size-independent local matrices to approximately 
but accurately represent the translation invariant ground 
states of all finite periodic systems. The performance 
of the algorithms for translational invariant systems ex- 
plored in this work is improved over the other algorithms 
by dropping the size dependence in the calculation of the 
wave function. The ground state wavefunctions of all fi- 
nite size systems can be obtained simply from the local 
matrices that are determined from an infinity lattice. 

To understand this, let us consider the following trans- 
lation invariant matrix product state in a ID bipartite 
system 

\9) = Tr(A[a 1 ]B[a 2 ]--'A[a L - 1 ]B[a L ])\a 1 "-a L ), (1) 

where \<ii) is the local basis state at site i and L is the 
lattice size. Given a, A[a] and B[a] are D x D matrices. 
The trace in Eq. |T]) ensures periodic boundary condition. 
In the ground state, if one determines the local matrices 
by minimizing the energy, then A and B should in general 
be L dependent. However, in many cases, such as in the 
valence bond solid state proposed by Affleck et al4, the 
local matrices can be L independent. In the following, 
we will show that one can always use matrix product 
states with L-indcpcndcnt A and B to approximately 
but accurately represent the ground state wavefunctions 
in finite periodic systems. 

The matrix product wavefunction of the ground state 
can be found by applying the projection operator 



2 



exp(—/3H) to an arbitrary initial state |\&o)- 

\V g ) = lim exp(-/3H)|* ). (2) 

This projection, which is equivalent to taking an imagi- 
nary time evolution, can be done iteratively. To do this, 
we divide /? into N steps with a small incremental time 
r = /3/N and decompose the projection operator at each 
step using the Trotter-Suzuki formula 

e -rH = e -rft„ e -rH eM „ + Q ^2^ (3) 

where 

-ffodd = J]] h,2i-l,2i, 
i 

Heven = h^iai+l- 

i 

At each iteration, the projection can be done in 
two successive steps, using the projection operators 
cxp(— rHodd) and exp(— rH even ), respectively. Since all 
the terms in H dd or H even commute with each other, 
these projections can be done by performing only local 
operations. For example, by applying exp(— rH dd) to 
Eq. ([1]), the matrix product wavefunction will become 

= Tr (M[a lt 0-2] ■ ■ ■ M[* L -i,<tl]) Wi ■ ■ ■ o- L ), (4) 
where M[a\,a 2 ] is a local matrix defined by 

M a p[ai,0-2] — Maurx^o-a 

= ]T A ai [a[]B~ <p [a' 2 ]{a 1 ,a 2 \e- Th ^\a[,cj l 2 ). (5) 

By singular value decomposition, one can decompose M 
to 

M a<Tlt0(72 = 2J U^ jXfV^^ , (6) 
l 

where Um and Vm are unitary matrices, A a/ is a semi- 
positive vector. From this we can rewrite M as a product 
of two matrices 

M[a 1 ,a 2 ]=A'[a 1 ]B'[a 2 l (7) 

KM\ = U^ l (\f I ) 1/ \ (8) 

^M = (Af) 1/2 ^ 2 . (9) 

is D\ x dD 2 matrix and B'[a 2 ] is a dD 2 x D\ 
matrix. Thus the wavefunction after the projection has 
the same matrix- product form as \^S). The only difference 
is that A and B in Eq. ([T]) are now replaced by A' and 
B', respectively. Given A and B, clearly the values of A' 
and B' are determined purely by the local Hamiltonian 
hi2 no matter how many sites in the system. 

In the limit r —> and /3 — Nt — > oo, the matrix 
product state obtained through above projection should 



approach the exact ground state. If there is no truncation 
to the matrix dimension, then the final site matrices A 
and B such obtained should be size independent. This 
means at least in the limit the bond dimension being 
infinity, the matrix product wavefunction of the ground 
state can be represented by the same local matrices, A 
and B, no matter how large the system size is. 

Of course, in practical calculation, one has to truncate 
the matrix dimension in order to carry out the projec- 
tion for sufficiently many times. Otherwise, the matrix 
dimension will blow up exponentially with the projec- 
tion steps. Upon truncation, the matrix product state 
with size-independent A and B will generally not be 
the rigorous ground state wavefunction. Nevertheless, it 
should still be a good approximation to the true ground 
state wavefunction. In particular, in the limit that the 
bond dimension approaches infinity, it should approach 
to the exact result. In the projection methods, the ac- 
curacy is limited only by the Trotter error controlled by 
t and the truncation error controlled by the matrix di- 
mension. 

The above argument can be readily extended to 
the tensor-network states in two or higher dimensions. 
Tensor-network wavefunctions, for example, can be de- 
termined using the entanglement mean-field projection 
approach introduced in Refs. @ and [l(|. Thus accurate 
ground state wave functions for all finite systems can be 
constructed by a few size independent local tensors. In 
higher dimensional systems, the lattice sizes that can be 
treated are generally very small due to the rapid growth 
of Hilbert space with the system size. The simplification 
of the problems to find a few tensors as suggested in this 
study is a big step towards solving the problems in two 
or higher dimensions. 

Besides the projection method, the matrix product 
state in one dimension can be also obtained using several 
other methods^. In an infinite lattice system, the local 
matrices obtained with different methods are equivalent. 
They can be gauged into a canonical formii by certain 
unitary transformations. 



II. DETERMINATION OF LOCAL MATRICES 

In this section, we introduce two approaches to eval- 
uate the local matrices. One is an entanglement mean 
field projection approach. This is an approach that was 
first proposed for evaluating a matrix product state in 
an infinite lattice via the imaginary time evolution in ID 
(Ref. [3) and 2D (Refs. d and OS). The other is the 
standard DMRG method with open boundary condition. 
This approach is applicable only in ID. 

For the discussion below, we will take the Heisenberg 
spin chain as an example to show how accurate a transla- 
tion invariant matrix product state with size independent 
local matrices can be. It is straightforward to extend the 
qualitative results to other quantum lattice models with 
short range interactions. The Hamiltonian of the Heisen- 
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berg model is defined by 

H = /lz,i+i, = Si ■ Si+i, (10) 

i 

where Si is the SU(2) spin operator. 

A. By projection 

Since the local matrices A and B are size independent, 
we can always use the methods that have been developed 
for studying a matrix product state in an infinite lattice 
to evaluate these matrices. In an infinite lattice, open 
boundary condition can be imposed without breaking the 
translation invariance of the matrix product state. In this 
case the two ends of the lattice are always disentangled 
and it is sufficient to perform just local transformation 
to convert a matrix product state into its canonical form. 
This can dramatically reduce the computational cost in 
the determination of local matrices. 

In this regards, a commonly adopted approach is to 
take all matrix elements of A and B as variational param- 
eters and determine them by minimizing the ground state 
energy. This, as discussed in Refs. 0-[l|, can be achieved 
by solving a generalized eigenvalue problem. Since there 
is no entanglement between the two ends of the lattice, 
the dimension of the boundary matrix reduces to 1 and 
the cost of this algorithm^ will scale as D 3 . 

Another kind of approach is to find the local matrices 
by performing the imaginary time evolution as defined 
by Eqs. (J2J7J) . These equations set a general framework 
for evaluating a matrix product state if no truncation 
is needed. However, in practical calculation, we have to 
truncate the matrix dimension at each step of projection. 
If we truncate the matrix dimension just using the singu- 
lar value decomposition of M-matrix defined by Eq. |(7J), 
the matrix product state will generally not converge to 
the true ground state. This is because M is purely a local 
matrix and the contribution from the environment ma- 
trix has not been considered in the basis truncatio n 10 ! 13 . 

The interplay between M and the environment matrix 
can be properly and accurately handled by performing a 
number of transformations to enable A and B to satisfy 



the following canonical conditions 
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where A Q and A& are semipositive diagonal matrices. The 
square of the diagonal matrix element of A a or is the 
eigenvalue of the reduced density matrix for the two semi- 
infinity blocks separated by AB or BA bond. It measures 



the entanglement between these two blocks and is equal 
to the possibility of the corresponding basis vector in the 
matrix product state. The truncation can be done by just 
keeping the largest D eigenvalues of A a ,b and the corre- 
sponding basis vectors. This can minimize the truncation 
error at each project step^, same as in the DMRG. The 
matrix product state such obtained is believed to be the 
most accurate one within the Trotter error. The cost 
of this algorithm also scales as D 3 . But this method is 
not applicable to a tensor-network state in two or higher 
dimensions. 

A more efficient and easy to implement approach, 
which will be used in the calculation below, is to evaluate 
the environment contribution by taking a mean field ap- 
proximation. This is a generalization of the poor-man's 
approach of second renormalization of tensor-network 
states introduced in Ref. [T^. The mean field parame- 
ters are the bond vectors that are introduced to approx- 
imately measure the entanglement between the matrix 
to be decomposed and rest of other matrices. No canoni- 
cal transformation needs to be done in this algorithm. At 
each step of projection, the truncation error is larger than 
that obtained by the canonical transformation. However, 
the truncation error is not accumulated in the process 
of iterations. The results obtained with this approach 
can reach the same accuracy as those obtained by the 
canonical transformation. 10 This approach converges fast 
and is highly reliable, provided that the short imaginary 
time evolution operator is sufficiently close to a unitary 
operator^ It can be used not just for studying matrix- 
product states in one dimension^, but also for studying 
tensor-network states in two or higher dimensions ^ 10 ' 13 . 
A detailed introduction to this method can be found from 
Ref. [121 ] for one-dimensional and Refs. for two- 

dimensional systems, respectively. 

To consider the renormalization effect of environment 
on M-matrix, let us redefine A and B as 

A[a] = T>]A a , (15) 
B[a] = r"[a]A b , (16) 

where A a and A& are positive diagonal matrices (also 
called bond vectors) defined on the bonds connecting two 
neighboring sites. \ a .b do not depend on a. In the canon- 
ical representation of the matrix product state, they are 
just the diagonal matrices defined in Eqs. (fT3"|) and (fl4"|) . 
In general, they can be considered as an approximation to 
the diagonal matrices in the canonical form. They mea- 
sure the entanglement between the left and right blocks 
connected by A a or Ah in an infinite system. The gauge 
degrees of freedom of the matrix product are also par- 
tially fixed by these bond vectors. 

If the basis states of the left and right blocks connected 
to the M-matrix defined by Eq. ([7]) are orthonormal, as in 
a canonical form, then the environmental contribution to 
M is proportional to Ah from left side of the environment 
and 1 from the right side of the environment. Thus the 
entanglement field on the bond for M is given by the 
singular values of the following environment modified M- 
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FIG. 1. The relative error of the ground state energy 
(E — Eo)/Eo as a function of the matrix dimension D for 
the spin half Heisenberg model with L = 30 and L = 100 
sites, respectively. Eq is the exact ground state energy. The 
results for other system sizes are similar. 



FIG. 2. The size dependence of the relative error of the ground 
state energy for the S=l/2 Heisenberg spin chain obtained 
with the matrix product state of bond dimension D = 50. 



matrix 

M[a 1 ,a 2 ] = XbM[a- 1 ,a 2 ] 

= X b T a [a 1 ]X a T b [a 2 ]X b . (17) 

In the mean-field calculation, however, the basis states 
on both the left and right blocks are not orthogonalized. 
In this case, A Qj b is just an approximate measure of the 
entanglement. Nevertheless, it still provides an efficient 
and good account of the renormalization effect of the 
environment. The reason for this is that in the projec- 
tion method, the truncation error is not accumulated and 
there is no need to minimize the truncation error at every 
step of iteration. 

By SVD, we can decompose M as a product of matrices 

M a ,p[<Jl,<T2] = Maa y ,Pa 3 = Uaa.l^lVl^a 2 , (18) 

where U and V are unitary matrices and A is a semi- 
positive diagonal matrix. If M[<7i,<72] is D x D matrix 
for given <j\ and a 2 , then U, A and V are all dD x dD 
matrices. To proceed the iteration, the matrix dimen- 
sion needs to be truncated. This can be done by keeping 
the largest D singular values of A and the corresponding 
vectors. After this one can rewrite the matrix-product 
wavefunction back to the form defined by Eq. ([T]) by up- 
dating r a ' b and A a with the formula 

XaM = Ay. 

Xb remains unchanged. This completes one step of imag- 
inary time evolution with H dd- Similarly one can carry 
out the imaginary time evolution with H even and update 



all site matrices and bond vectors. The converged ma- 
trices from the iterated projections are used to construct 
ground state wavefunctions for all system sizes. 

The error comes from two sources: one is the trun- 
cation error and the other is the Trotter error arisen 
from the Trotter-Suzuki decomposition. The Trotter er- 
ror does not depend on the system size and can be re- 
duced by using a smaller r or using higher order Trotter- 
Suzuki decomposition formula. The truncation error is 
controlled by the matrix dimension D. It can be reduced 
by increasing the matrix dimension. 

Fig. [T] shows the relative error of the ground state en- 
ergy as a function of the matrix dimension D for the 
spin- 1/2 Heisenberg model on finite lattice systems ob- 
tained from the matrix product wavefunction with size- 
independent A and B. We only show the results of rel- 
atively small system sizes, as the success of this method 
in the large system size limit is well documented. In the 
regime of small D, the truncation error dominates, the 
error drops quickly with increasing D. The error stops 
dropping when D becomes larger than certain critical 
bond dimension D c , beyond which the truncation error 
becomes smaller than the Trotter error. The value of D c 
depends on the small time interval r and the order of 
Trotter-Suzuki decomposition formula used in the imag- 
inary time evolution. But it does not depend much on 
the system size. This is probably because the same local 
matrices are used in the calculation of the ground state 
energy no matter how large the system size is. In ob- 
taining the results in Fig. [T] we have used the first order 
Trotter-Suzuki decomposition formula. In this case, the 
Trotter error scales approximately as r 2 . This is consis- 
tent with the numerical results shown in Fig. [T] 

Fig. [5] shows the lattice size dependence of the relative 
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FIG. 3. The relative error of the ground state energy, (E — 
Eo)/Eo, as a function of matrix dimension D for the spin one 
Heisenberg model. Eo is the exact ground state energy. 

error of the ground state of the S=l/2 Heisenberg model 
with D = 50. The values of Eq are from the exact diag- 
onalization results published in Ref. [15j . For the given 
t, D — 50 are larger than the critical D c , the errors are 
controlled by the Trotter error and almost do not depend 
on the system size. 

The matrix dimension dependence of the relative er- 
rors of the ground state energy behaves similarly for the 
S=l Heisenberg spin chain, as shown in Fig. |3l Unlike 
the S=l/2 Heisenberg model whose excitation spectrum 
is critical, the S=l Heisenberg model has a finite excita- 
tion gap and its ground state energy converges quickly 
with increasing D. Again the error is dominated by the 
truncation error in the small D regime and by the Trotter 
error when D is larger than a r-dependent critical D c . 

The above results show unambiguously that the ma- 
trix product states constructed with size-independent lo- 
cal matrices provide accurate and unified description of 
ground states irrespective of the system size. The accu- 
racy of the wavefunction can be systematically and re- 
liably improved by increasing the matrix dimension D 
to reduce the truncation error and by using a smaller r 
or higher order Trotter-Suzuki decomposition formula to 
reduce the Trotter error. 



B. By the DMRG 

In a translation invariant system, one can also use the 
DMRG method for open boundary systems to find the 
local matrices A and B defined by Eq. ([IJ . To understand 
this, let us consider the ground state generated by the 
DMRG with open boundary condition 

I*) = E ^2^( s o^i^2,e 2 )\s a 1 a 2 e 2 ), (19) 



where we have used 1 and 2 to denote the coordinates 
of two middle sites of the chain, s ) and e 2 ) are the 
basis states of the system and environment blocks, re- 
spectively. By singular value decomposition, we can de- 
compose Vt^so, oi, o 2 , e 2 ) as 

= E< ) ^W 1) CW (20) 

I 

In the above expression, the singular value A^ 1 ' is noth- 
ing but the square root of the eigenvalue of the reduced 
density matrix. U and V are the unitary (or isometric af- 
ter truncation) basis transformation matrices. The basis 
states of system and environment blocks are transformed 
according to the following equations 

\ g i)= E UjV^Msi-iVi), (21) 
\ e *)= E y^e i+1 Wi+i]Wi+ie i+ i}. (22) 

e; + iCT i+ i 

Substituting Eq. ([20} into Eq. ([i"9]>. we obtain 

l*) = E E O^WN2Ml s o<no 2 e 2 ). (23) 

By using the basis transformation formula (l2~Tj) and (|22|l , 
we can further express this wavefunction as 

i*>= E E ^Ij^Smaw 

<T0tTiCT2tT3 S_iSQe2e3,i 

Vf% [o 2 ]K ( 2 3) e3 [03] |S-lO OxO 2 O3e3). (24) 

If the system is reflection symmetric with respect to the 
middle bond of the chain, and should be the 
Hermitian conjugate of and U(°\ respectively. In 
the DMRG iteration, the two local unitary matrices on 
the system block besides the middle bond, and U^ x \ 
will converge alternatively in a bipartite lattice model 
when the system size becomes sufficiently largei£. If 
we use these two converged isometric matrices to form 
a translation invariant matrix product state, i.e. 

A[a] = U^[<7) ) (25) 
B[a] = uW[a], (26) 

we expect that it will be a good approximation to the 
ground state of the Hamiltonian with periodic boundary 
condition. This is because in an infinite system, the wave- 
function in the middle of the chain should be the same 
no matter what kind of boundary condition is used. 

The local matrices such obtained are Trotter-Suzuki 
decomposition error free. Fig. 3] compares the relative er- 
ror of the ground state energy for the S =1/2 Heisenberg 
model obtained by the DMRG method to that obtained 
by the entanglement mean-field projection method. In 
the entanglement mean-field projection calculation, we 
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FIG. 4. Comparison of the relative error of the ground state 
energy of the spin-half Heisenberg model as a function of D 
obtained using the matrix product states whose local matrices 
are determined by the DMRG with those determined by the 
poor-man's projection. 



have chosen a small r so that the Trotter error is much 
smaller than the truncation error. We find that the er- 
rors obtained with these two methods are of the same 
order. 

To construct the translation invariant matrix product 
states using the DMRG is generally less efficient than 
the entanglement mean-field projection method. But the 
DMRG calculation can be combined with the conven- 
tional DMRG study with open boundary condition. No 
extra cost is needed in order to calculate the local matri- 
ces used by the matrix product state. 

In the above discussion, we have assumed that the iso- 
metric matrices J7(°) and f/W would converge for suffi- 
ciently large chain. This is literally correct. However, if 
there are degeneracy or nearly degeneracy within com- 
puter machine error in the singular values A in Eq. ([20]) . 



{/-matrices may not be uniquely fixed since the ground 
state is unchanged by swapping any pair of these degener- 
ate states. Consequently the column index of B[<r] may 
not perfectly match the row index of A[a] (there is no 
mismatch between the row index of B[a\ and the column 
index of A[<j]). In this case, the local matrices defined by 
Eqs. (|25H26|) may not be a good description of the ground 
state. 



III. SUMMARY 

We have shown that the local matrices (tensors) ob- 
tained from an infinite lattice can be also used to ac- 
curately represent the matrix product states (or tensor- 
network states) in finite translation invariant lattice sys- 
tems. This provides an efficient way to determine the 
matrix product states or tensor-network states in one- or 
higher-dimensional periodic systems. 

The size- independent local matrices can be determined 
by the approximate entanglement projection with a poor- 
man's treatment to the environment lattice. This method 
can be applied not only in ID, but also in two- or 
higher dimensions. The accuracy of the wavefunction 
such obtained is controlled by both the truncation and 
the Trotter-Suzuki decomposition errors. For small D, 
the truncation error dominates. This error, however, 
can be reduced below the Trotter- Suzuki decomposition 
error simply by increasing the bond dimension D. In 
that case, the accuracy of the matrix product wavefunc- 
tion is purely determined by the Trotter-Suzuki decom- 
position error. This error can be reduced by using a 
small Trotter-Suzuki parameter r or using a higher or- 
der Trotter-Suzuki decomposition formula. One can also 
find the size-independent local matrices using the DMRG 
with open boundary conditions. The matrix product 
states such obtained do not have the Trotter-Suzuki de- 
composition error. 

We acknowledge the support of NSF-China and the 
National Program for Basic Research of the Ministry of 
Science and Technology of China. 
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